function quadrupole
  close all
  Omega = 1.1; omega = 10; beta = 20;
  tspan = [0:0.2:100];
  initial_w = [0.005,0.005,0,0,0,0.01];
  [times,sols] = ode45(@(t,w) eom(t,w,Omega,omega,beta),tspan,initial_w);  
  animate_quadrupole(times,sols); 
%   solution to (c)
  for i=1:101
      Omega(i) = 0.6 + 1.1*(i-1)/100;
     [times,sols] = ode45(@(t,w) eom(t,w,Omega(i),omega,beta),tspan,initial_w);
     rval(i) = sqrt(sols(end,1)^2+sols(end,2)^2);
  end
  figure
  plot(Omega,rval,'LineWidth',2,'Color','r')
  ylim([0,1]);
  set(gca,'FontSize',14);
  xlabel('\Omega');
  ylabel('dist from quadrupole axis (m)');
  title('Solution to (c)');  
end
function dwdt = eom(t,w,Omega,omega,beta)
  position = w(1:3); velocity = w(4:6);
  accel = Omega^2*(1+beta*cos(omega*t))*[position(1);position(2);0];
  dwdt = [velocity;accel];
end
function animate_quadrupole(times,sols)

    [XX,YY,ZZ] = cylinder(0.1);
    XX1 = XX - 0.15;
    YY1 = YY - 0.15;
    XX2 = XX + 0.15;
    YY2 = YY + 0.15;
    figure1 = figure;
    for i=1:length(times)
        clf
     axes1 = axes('Parent',figure1,'DataAspectRatio',[1 1 1]);

     view(axes1,[-79.5 22]);
     grid(axes1,'on');
     xlim([-0.6 0.6]);
     ylim([-0.6 0.6]);
     zlim([0 1.2]);
     hold on
     surf(XX1,YY1,ZZ)
     surf(XX2,YY2,ZZ)
     surf(XX1,YY2,ZZ)
     surf(XX2,YY1,ZZ)      
     plot3(sols(1:i,1),sols(1:i,2),sols(1:i,3),...
         'LineWidth',2,'Color',[0 1 0]);
     plot3(sols(i,1),sols(i,2),sols(i,3),...
         'ro','MarkerSize',11,'MarkerFaceColor','r');
     pause(0.01)
    end
end